
/* ----------------------------------------------------------------------
 * G-Nut - GNSS software development library
 * 
  (c) 2018 G-Nut Software s.r.o. (software@gnutsoftware.com)
  This file is part of the G-Nut C++ library.
 
-*/

#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <algorithm>

#include "newmat/newmatap.h"

#include "gutils/gsysconv.h"
#include "gmodels/gtide2010.h"
#include "gutils/gmatrixconv.h"
#include "gmodels/ginterp.h"

using namespace std;

namespace gnut {
// constructor
// ----------
t_gtide2010::t_gtide2010(t_glog* l) : t_gtide(l)
{
   _D.ReSize(6);
   _DD.ReSize(6);   
}


// destructor
// ----------
t_gtide2010::~t_gtide2010()
{}

// WARNING!!! Still IERS 2003. 2010 needs to be implemented
// solid earth tides
// ----------
t_gtriple t_gtide2010::tide_searth(const t_gtime& epoch, t_gtriple& crd)
{
#ifdef BMUTEX   
   boost::mutex::scoped_lock lock(_mutex);
#endif 
   _mutex.lock();

  t_gtriple dxyz(0.0,0.0,0.0);

  ColumnVector xSun;
  ColumnVector xMoon;  
  double       rSun;
  double       rMoon;

  double Mjd = epoch.dmjd();

  xSun = _planEph.sunPos(Mjd).crd_cvect();
//  cout << "XSun " << xSun << endl;
  rSun = sqrt(DotProduct(xSun,xSun));
  xSun /= rSun;
  xMoon = _planEph.moonPos(Mjd).crd_cvect();
  rMoon = sqrt(DotProduct(xMoon,xMoon));
  xMoon /= rMoon;
   
  ColumnVector xyz = crd.crd_cvect();   
  double       rRec    = sqrt(DotProduct(xyz, xyz));
  ColumnVector xyzUnit = xyz / rRec;
  t_gtriple ell;
  xyz2ell(crd, ell, false);

  // Love's, Shida's Numbers
  // --------------
  const double    h0 =  0.6078;   // Love number h0 of degree 2
  const double    h2 = -0.0006;   // Love number h2 of degree 2
  const double    h3 =  0.292;    // Love number of degree 3   
  const double    l0 =  0.0847;   // Shida number l0 of degree 2
  const double    l2 =  0.0002;   // Shida number l0 of degree 2
  const double    l3 =  0.015;    // Shide number of degree 3 
  
  double h = h0 + h2*(3*sin(ell[0])*sin(ell[0]) - 1)/2;
  double l = l0 + l2*(3*sin(ell[0])*sin(ell[0]) - 1)/2;
   
  const double    MS2E = 332946.0;        // Ratio of mass Sun to Earth
  const double    MM2E = 0.01230002;      // Ratio of mass Moon to Earth
      
  // Tidal Displacement
  // ------------------
  double scSunRec  = DotProduct(xSun,  xyzUnit);
  double scMoonRec = DotProduct(xMoon, xyzUnit);
   
  // Computing temporary vectors
  t_gtriple tmpSun, tmpMoon;  
  for (int i = 0 ; i < 3; i++) {
     tmpSun[i]  = xSun[i]  - scSunRec *xyzUnit[i];
     tmpMoon[i] = xMoon[i] - scMoonRec*xyzUnit[i];
  }
      
  for(int i = 0; i < 3; i++){
     dxyz[i] = MS2E*(pow(rRec,4)/pow(rSun,3)) *(h*xyzUnit[i]*(3.0/2.0*scSunRec *scSunRec  - 0.5) + 3.0*l*scSunRec *tmpSun[i]) +
               MM2E*(pow(rRec,4)/pow(rMoon,3))*(h*xyzUnit[i]*(3.0/2.0*scMoonRec*scMoonRec - 0.5) + 3.0*l*scMoonRec*tmpMoon[i]);
          
     dxyz[i] += MM2E * (pow(rRec,5)/pow(rMoon,4))*
                                    (h3 * xyzUnit[i] * (5.0/2.0*scMoonRec*scMoonRec*scMoonRec - 3.0/2.0*scMoonRec) +
                                    l3 * (15.0/2.0*scMoonRec*scMoonRec - 3.0/2.0) * tmpMoon[i]);
  }
     
//cout << "return tide: " << dxyz[0] << " " << dxyz[1] << " " << dxyz[2] << endl;
  _mutex.unlock(); 
   return dxyz;
}


// pole tides
// ----------
t_gtriple t_gtide2010::tide_pole()
{
#ifdef BMUTEX   
   boost::mutex::scoped_lock lock(_mutex);
#endif 
   _mutex.lock();

  t_gtriple dxyz(0.0,0.0,0.0);

  _mutex.unlock();
  return dxyz;
}


// ---------------------------------------------------------------------
// <0, 2pi)
// ---------------------------------------------------------------------
double null2pi_degree(double angle)
{
    while (angle > 360.0)
        angle -= 360.0;
    while (angle < 0)
        angle += 360.0;
    return angle;
}

double null2pi_rad(double angle)
{
    while (angle > 2*G_PI)
        angle -=  2*G_PI;
    while (angle < 0)
        angle +=  2*G_PI;
    return angle;
}

// ---------------------------------------------------------------------
// Astronomical arguments and fundamental frequencies
// ---------------------------------------------------------------------
Matrix DD_arguments(double mjd, double delta)
 {
    Matrix Kij(5,5), tM(5,1);                                        // general
    Matrix Delauney(6,1), l(1,5), ls(1,5), F(1,5), D(1,5), Om(1,5), fund_freqDel(5,1); // Delauney
    Matrix dood(6,3), Doodson0(6,5), Doodson(6,5), Doodson_vec(6,1), Doodson_vecf(12,1), s(1,5), h(1,5), p(1,5), Ns(1,5), ps(1,5), T(1,5), fund_freqDood(6,1);  // Doodson
    double t0, t1, t2, t3, t4;

	// nulling
	dood = 0.0;
	Doodson0 = 0.0;
	Delauney = 0.0;
	Doodson_vec=0.0;

    t0=1.0;
    t1=(mjd - 51544.5 + delta/86400.0)/36525.0; // 51544.5 vs. 51545.0.... double T = (dmjd - 51545.0 + delta/86400.0)/36525.0; 
    t2=pow (t1,2);
    t3=pow (t1,3);
    t4=pow (t1,4);
    tM << t0 << t1 << t2 << t3 << t4;
    
    // IERS 2010, p.67 ... l, l', F, D, Omega (F1 ... F5)
    Kij << 134.96340251         << 357.52910918        << 93.27209062         << 297.85019547         <<  125.04455501
		<< 1717915923.2178/3600 << 129596581.0481/3600 << 1739527262.8478/3600<< 1602961601.2090/3600 << -6962890.5431/3600 // /3600;
		<< 31.8792/3600         <<-0.5532 /3600        <<-12.7512/3600        << -6.3706/3600         <<  7.4722/3600 // /3600;
		<<  0.051635/3600       << 0.000136/3600       <<-0.001037/3600       <<  0.006593/3600       <<  0.007702/3600 // /3600;
	    << -0.00024470/3600     <<-0.00001149/3600     << 0.00000417/3600     << -0.00003169/3600     << -0.00005939/3600; // /3600;

    fund_freqDel << Kij(2,1) << Kij(2,2) << Kij(2,3) << Kij(2,4) << Kij(2,5);// in [degree/Julcent]

	// Delauney arguments, IERS 2010 p.48
	l << Kij(1,1) << Kij(2,1) << Kij(3,1) << Kij(4,1) << Kij(5,1);
	ls<< Kij(1,2) << Kij(2,2) << Kij(3,2) << Kij(4,2) << Kij(5,2);
	F << Kij(1,3) << Kij(2,3) << Kij(3,3) << Kij(4,3) << Kij(5,3);
	D << Kij(1,4) << Kij(2,4) << Kij(3,4) << Kij(4,4) << Kij(5,4);
	Om<< Kij(1,5) << Kij(2,5) << Kij(3,5) << Kij(4,5) << Kij(5,5);
  
	//  Doodson arguments, Dehant (2015) et al., p.178
	s  = F + Om;          // mean lunar node
	h  = F + Om - D;      // angle of the mean tropic year
	p  = -l + F + Om;     // mean lunar perigee
	Ns = -Om;             // angle of the mean lunar node
	ps = -ls + F + Om -D; // angle of perihelion - ps/p'
	T = - s; // -s [+ gmst + pi] ... added below
	
	// Delauney
	Delauney << DotProduct(T,tM) << DotProduct(l,tM) << DotProduct(ls,tM) << DotProduct(F,tM) << DotProduct(D,tM) << DotProduct(Om,tM);
	Delauney=Delauney*D2R;
	
	// Doodson0 ... T is not ready here!
	Doodson0.Row(1) = T;
	Doodson0.Row(2) = s;
	Doodson0.Row(3) = h;
	Doodson0.Row(4) = p;
	Doodson0.Row(5) = Ns;
	Doodson0.Row(6) = ps;

	// Doodson fundamental frequencies
	fund_freqDood = Doodson0.Column(2)/36525.0/24.0; // in [degree/hour] in mean Solar day
	fund_freqDood.Row(1) = 15.0 - fund_freqDood.Row(2) + fund_freqDood.Row(3);

	// Doodson ... final modification
	Doodson = Doodson0;
	// convert to <0,360)
	Doodson(1,1)= null2pi_degree(Doodson(1,1));
	Doodson(2,1)= null2pi_degree(Doodson(2,1));
	Doodson(3,1)= null2pi_degree(Doodson(3,1));
	Doodson(4,1)= null2pi_degree(Doodson(4,1));
	Doodson(5,1)= null2pi_degree(Doodson(5,1));
	Doodson(6,1)= null2pi_degree(Doodson(6,1));
	
	// make a vector
	Doodson_vec = Doodson * tM;
	
	// get GMST
	t_gephplan epoGMST;
	double gmst_i = R2D*epoGMST.gmst(mjd);
	
	// update T for [+ gmst + pi]!
	Doodson_vec(1,1) = Doodson_vec(1,1) + gmst_i + 180.0;
	
	// convert to <0,360) and put to Doodson_vecf
	Doodson_vecf(1,1)= D2R*null2pi_degree(Doodson_vec(1,1)); // zde se diky starsi funkci GMST v gephplan.h lisime na 8 des. miste (nema vliv na vysledek)
	Doodson_vecf(2,1)= D2R*null2pi_degree(Doodson_vec(2,1));
	Doodson_vecf(3,1)= D2R*null2pi_degree(Doodson_vec(3,1));
	Doodson_vecf(4,1)= D2R*null2pi_degree(Doodson_vec(4,1));
	Doodson_vecf(5,1)= D2R*null2pi_degree(Doodson_vec(5,1));
	Doodson_vecf(6,1)= D2R*null2pi_degree(Doodson_vec(6,1));
	
	// add frequencies
	Doodson_vecf(7,1)=fund_freqDood(1,1);
	Doodson_vecf(8,1)=fund_freqDood(2,1);
	Doodson_vecf(9,1)=fund_freqDood(3,1);
	Doodson_vecf(10,1)=fund_freqDood(4,1);
	Doodson_vecf(11,1)=fund_freqDood(5,1);
	Doodson_vecf(12,1)=fund_freqDood(6,1);
	
    return Doodson_vecf;
  }

Matrix dnn_table()
 {
    Matrix dnntab(342,7); 
    dnntab // first 6. cols = Doodson number, 7. col = tide amplitude from TGP
    <<     2 <<    0  <<   0  <<   0  <<   0  <<   0 <<  0.632208
    <<     2 <<    2  <<  -2  <<   0  <<   0  <<   0 <<  0.294107
    <<     2 <<   -1  <<   0  <<   1  <<   0  <<   0 <<  0.121046
    <<     2 <<    2  <<   0  <<   0  <<   0  <<   0 <<  0.079915
    <<     2 <<    2  <<   0  <<   0  <<   1  <<   0 <<  0.023818
    <<     2 <<    0  <<   0  <<   0  <<  -1  <<   0 << -0.023589
    <<     2 <<   -1  <<   2  <<  -1  <<   0  <<   0 <<  0.022994
    <<     2 <<   -2  <<   2  <<   0  <<   0  <<   0 <<  0.019333
    <<     2 <<    1  <<   0  <<  -1  <<   0  <<   0 << -0.017871
    <<     2 <<    2  <<  -3  <<   0  <<   0  <<   1 <<  0.017192
    <<     2 <<   -2  <<   0  <<   2  <<   0  <<   0 <<  0.016018
    <<     2 <<   -3  <<   2  <<   1  <<   0  <<   0 <<  0.004671
    <<     2 <<    1  <<  -2  <<   1  <<   0  <<   0 << -0.004662
    <<     2 <<   -1  <<   0  <<   1  <<  -1  <<   0 << -0.004519
    <<     2 <<    3  <<   0  <<  -1  <<   0  <<   0 <<  0.004470
    <<     2 <<    1  <<   0  <<   1  <<   0  <<   0 <<  0.004467
    <<     2 <<    2  <<   0  <<   0  <<   2  <<   0 <<  0.002589
    <<     2 <<    2  <<  -1  <<   0  <<   0  <<  -1 << -0.002455
    <<     2 <<    0  <<  -1  <<   0  <<   0  <<   1 << -0.002172
    <<     2 <<    1  <<   0  <<   1  <<   1  <<   0 <<  0.001972
    <<     2 <<    3  <<   0  <<  -1  <<   1  <<   0 <<  0.001947
    <<     2 <<    0  <<   1  <<   0  <<   0  <<  -1 <<  0.001914
    <<     2 <<    0  <<  -2  <<   2  <<   0  <<   0 << -0.001898
    <<     2 <<   -3  <<   0  <<   3  <<   0  <<   0 <<  0.001802
    <<     2 <<   -2  <<   3  <<   0  <<   0  <<  -1 <<  0.001304
    <<     2 <<    4  <<   0  <<   0  <<   0  <<   0 <<  0.001170
    <<     2 <<   -1  <<   1  <<   1  <<   0  <<  -1 <<  0.001130
    <<     2 <<   -1  <<   3  <<  -1  <<   0  <<  -1 <<  0.001061
    <<     2 <<    2  <<   0  <<   0  <<  -1  <<   0 << -0.001022
    <<     2 <<   -1  <<  -1  <<   1  <<   0  <<   1 << -0.001017
    <<     2 <<    4  <<   0  <<   0  <<   1  <<   0 <<  0.001014
    <<     2 <<   -3  <<   4  <<  -1  <<   0  <<   0 <<  0.000901
    <<     2 <<   -1  <<   2  <<  -1  <<  -1  <<   0 << -0.000857
    <<     2 <<    3  <<  -2  <<   1  <<   0  <<   0 <<  0.000855
    <<     2 <<    1  <<   2  <<  -1  <<   0  <<   0 <<  0.000855
    <<     2 <<   -4  <<   2  <<   2  <<   0  <<   0 <<  0.000772
    <<     2 <<    4  <<  -2  <<   0  <<   0  <<   0 <<  0.000741
    <<     2 <<    0  <<   2  <<   0  <<   0  <<   0 <<  0.000741
    <<     2 <<   -2  <<   2  <<   0  <<  -1  <<   0 << -0.000721
    <<     2 <<    2  <<  -4  <<   0  <<   0  <<   2 <<  0.000698
    <<     2 <<    2  <<  -2  <<   0  <<  -1  <<   0 <<  0.000658
    <<     2 <<    1  <<   0  <<  -1  <<  -1  <<   0 <<  0.000654
    <<     2 <<   -1  <<   1  <<   0  <<   0  <<   0 << -0.000653
    <<     2 <<    2  <<  -1  <<   0  <<   0  <<   1 <<  0.000633
    <<     2 <<    2  <<   1  <<   0  <<   0  <<  -1 <<  0.000626
    <<     2 <<   -2  <<   0  <<   2  <<  -1  <<   0 << -0.000598
    <<     2 <<   -2  <<   4  <<  -2  <<   0  <<   0 <<  0.000590
    <<     2 <<    2  <<   2  <<   0  <<   0  <<   0 <<  0.000544
    <<     2 <<   -4  <<   4  <<   0  <<   0  <<   0 <<  0.000479
    <<     2 <<   -1  <<   0  <<  -1  <<  -2  <<   0 << -0.000464
    <<     2 <<    1  <<   2  <<  -1  <<   1  <<   0 <<  0.000413
    <<     2 <<   -1  <<  -2  <<   3  <<   0  <<   0 << -0.000390
    <<     2 <<    3  <<  -2  <<   1  <<   1  <<   0 <<  0.000373
    <<     2 <<    4  <<   0  <<  -2  <<   0  <<   0 <<  0.000366
    <<     2 <<    0  <<   0  <<   2  <<   0  <<   0 <<  0.000366
    <<     2 <<    0  <<   2  <<  -2  <<   0  <<   0 << -0.000360
    <<     2 <<    0  <<   2  <<   0  <<   1  <<   0 << -0.000355
    <<     2 <<   -3  <<   3  <<   1  <<   0  <<  -1 <<  0.000354
    <<     2 <<    0  <<   0  <<   0  <<  -2  <<   0 <<  0.000329
    <<     2 <<    4  <<   0  <<   0  <<   2  <<   0 <<  0.000328
    <<     2 <<    4  <<  -2  <<   0  <<   1  <<   0 <<  0.000319
    <<     2 <<    0  <<   0  <<   0  <<   0  <<   2 <<  0.000302
    <<     2 <<    1  <<   0  <<   1  <<   2  <<   0 <<  0.000279
    <<     2 <<    0  <<  -2  <<   0  <<  -2  <<   0 << -0.000274
    <<     2 <<   -2  <<   1  <<   0  <<   0  <<   1 << -0.000272
    <<     2 <<   -2  <<   1  <<   2  <<   0  <<  -1 <<  0.000248
    <<     2 <<   -1  <<   1  <<  -1  <<   0  <<   1 << -0.000225
    <<     2 <<    5  <<   0  <<  -1  <<   0  <<   0 <<  0.000224
    <<     2 <<    1  <<  -3  <<   1  <<   0  <<   1 << -0.000223
    <<     2 <<   -2  <<  -1  <<   2  <<   0  <<   1 << -0.000216
    <<     2 <<    3  <<   0  <<  -1  <<   2  <<   0 <<  0.000211
    <<     2 <<    1  <<  -2  <<   1  <<  -1  <<   0 <<  0.000209
    <<     2 <<    5  <<   0  <<  -1  <<   1  <<   0 <<  0.000194
    <<     2 <<   -4  <<   0  <<   4  <<   0  <<   0 <<  0.000185
    <<     2 <<   -3  <<   2  <<   1  <<  -1  <<   0 << -0.000174
    <<     2 <<   -2  <<   1  <<   1  <<   0  <<   0 << -0.000171
    <<     2 <<    4  <<   0  <<  -2  <<   1  <<   0 <<  0.000159
    <<     2 <<    0  <<   0  <<   2  <<   1  <<   0 <<  0.000131
    <<     2 <<   -5  <<   4  <<   1  <<   0  <<   0 <<  0.000127
    <<     2 <<    0  <<   2  <<   0  <<   2  <<   0 <<  0.000120
    <<     2 <<   -1  <<   2  <<   1  <<   0  <<   0 <<  0.000118
    <<     2 <<    5  <<  -2  <<  -1  <<   0  <<   0 <<  0.000117
    <<     2 <<    1  <<  -1  <<   0  <<   0  <<   0 <<  0.000108
    <<     2 <<    2  <<  -2  <<   0  <<   0  <<   2 <<  0.000107
    <<     2 <<   -5  <<   2  <<   3  <<   0  <<   0 <<  0.000105
    <<     2 <<   -1  <<  -2  <<   1  <<  -2  <<   0 << -0.000102
    <<     2 <<   -3  <<   5  <<  -1  <<   0  <<  -1 <<  0.000102
    <<     2 <<   -1  <<   0  <<   0  <<   0  <<   1 <<  0.000099
    <<     2 <<   -2  <<   0  <<   0  <<  -2  <<   0 << -0.000096
    <<     2 <<    0  <<  -1  <<   1  <<   0  <<   0 <<  0.000095
    <<     2 <<   -3  <<   1  <<   1  <<   0  <<   1 << -0.000089
    <<     2 <<    3  <<   0  <<  -1  <<  -1  <<   0 << -0.000085
    <<     2 <<    1  <<   0  <<   1  <<  -1  <<   0 << -0.000084
    <<     2 <<   -1  <<   2  <<   1  <<   1  <<   0 << -0.000081
    <<     2 <<    0  <<  -3  <<   2  <<   0  <<   1 << -0.000077
    <<     2 <<    1  <<  -1  <<  -1  <<   0  <<   1 << -0.000072
    <<     2 <<   -3  <<   0  <<   3  <<  -1  <<   0 << -0.000067
    <<     2 <<    0  <<  -2  <<   2  <<  -1  <<   0 <<  0.000066
    <<     2 <<   -4  <<   3  <<   2  <<   0  <<  -1 <<  0.000064
    <<     2 <<   -1  <<   0  <<   1  <<  -2  <<   0 <<  0.000063
    <<     2 <<    5  <<   0  <<  -1  <<   2  <<   0 <<  0.000063
    <<     2 <<   -4  <<   5  <<   0  <<   0  <<  -1 <<  0.000063
    <<     2 <<   -2  <<   4  <<   0  <<   0  <<  -2 <<  0.000062
    <<     2 <<   -1  <<   0  <<   1  <<   0  <<   2 <<  0.000062
    <<     2 <<   -2  <<  -2  <<   4  <<   0  <<   0 << -0.000060
    <<     2 <<    3  <<  -2  <<  -1  <<  -1  <<   0 <<  0.000056
    <<     2 <<   -2  <<   5  <<  -2  <<   0  <<  -1 <<  0.000053
    <<     2 <<    0  <<  -1  <<   0  <<  -1  <<   1 <<  0.000051
    <<     2 <<    5  <<  -2  <<  -1  <<   1  <<   0 <<  0.000050
    <<     1 <<    1  <<   0  <<   0  <<   0  <<   0 <<  0.368645
    <<     1 <<   -1  <<   0  <<   0  <<   0  <<   0 << -0.262232
    <<     1 <<    1  <<  -2  <<   0  <<   0  <<   0 << -0.121995
    <<     1 <<   -2  <<   0  <<   1  <<   0  <<   0 << -0.050208
    <<     1 <<    1  <<   0  <<   0  <<   1  <<   0 <<  0.050031
    <<     1 <<   -1  <<   0  <<   0  <<  -1  <<   0 << -0.049470
    <<     1 <<    2  <<   0  <<  -1  <<   0  <<   0 <<  0.020620
    <<     1 <<    0  <<   0  <<   1  <<   0  <<   0 <<  0.020613
    <<     1 <<    3  <<   0  <<   0  <<   0  <<   0 <<  0.011279
    <<     1 <<   -2  <<   2  <<  -1  <<   0  <<   0 << -0.009530
    <<     1 <<   -2  <<   0  <<   1  <<  -1  <<   0 << -0.009469
    <<     1 <<   -3  <<   2  <<   0  <<   0  <<   0 << -0.008012
    <<     1 <<    0  <<   0  <<  -1  <<   0  <<   0 <<  0.007414
    <<     1 <<    1  <<   0  <<   0  <<  -1  <<   0 << -0.007300
    <<     1 <<    3  <<   0  <<   0  <<   1  <<   0 <<  0.007227
    <<     1 <<    1  <<  -3  <<   0  <<   0  <<   1 << -0.007131
    <<     1 <<   -3  <<   0  <<   2  <<   0  <<   0 << -0.006644
    <<     1 <<    1  <<   2  <<   0  <<   0  <<   0 <<  0.005249
    <<     1 <<    0  <<   0  <<   1  <<   1  <<   0 <<  0.004137
    <<     1 <<    2  <<   0  <<  -1  <<   1  <<   0 <<  0.004087
    <<     1 <<    0  <<   2  <<  -1  <<   0  <<   0 <<  0.003944
    <<     1 <<    2  <<  -2  <<   1  <<   0  <<   0 <<  0.003943
    <<     1 <<    3  <<  -2  <<   0  <<   0  <<   0 <<  0.003420
    <<     1 <<   -1  <<   2  <<   0  <<   0  <<   0 <<  0.003418
    <<     1 <<    1  <<   1  <<   0  <<   0  <<  -1 <<  0.002885
    <<     1 <<    1  <<  -1  <<   0  <<   0  <<   1 <<  0.002884
    <<     1 <<    4  <<   0  <<  -1  <<   0  <<   0 <<  0.002160
    <<     1 <<   -4  <<   2  <<   1  <<   0  <<   0 << -0.001936
    <<     1 <<    0  <<  -2  <<   1  <<   0  <<   0 <<  0.001934
    <<     1 <<   -2  <<   2  <<  -1  <<  -1  <<   0 << -0.001798
    <<     1 <<    3  <<   0  <<  -2  <<   0  <<   0 <<  0.001690
    <<     1 <<   -1  <<   0  <<   2  <<   0  <<   0 <<  0.001689
    <<     1 <<   -1  <<   0  <<   0  <<  -2  <<   0 <<  0.001516
    <<     1 <<    3  <<   0  <<   0  <<   2  <<   0 <<  0.001514
    <<     1 <<   -3  <<   2  <<   0  <<  -1  <<   0 << -0.001511
    <<     1 <<    4  <<   0  <<  -1  <<   1  <<   0 <<  0.001383
    <<     1 <<    0  <<   0  <<  -1  <<  -1  <<   0 <<  0.001372
    <<     1 <<    1  <<  -2  <<   0  <<  -1  <<   0 <<  0.001371
    <<     1 <<   -3  <<   0  <<   2  <<  -1  <<   0 << -0.001253
    <<     1 <<    1  <<   0  <<   0  <<   2  <<   0 << -0.001075
    <<     1 <<    1  <<  -1  <<   0  <<   0  <<  -1 <<  0.001020
    <<     1 <<   -1  <<  -1  <<   0  <<   0  <<   1 <<  0.000901
    <<     1 <<    0  <<   2  <<  -1  <<   1  <<   0 <<  0.000865
    <<     1 <<   -1  <<   1  <<   0  <<   0  <<  -1 << -0.000794
    <<     1 <<   -1  <<  -2  <<   2  <<   0  <<   0 <<  0.000788
    <<     1 <<    2  <<  -2  <<   1  <<   1  <<   0 <<  0.000782
    <<     1 <<   -4  <<   0  <<   3  <<   0  <<   0 << -0.000747
    <<     1 <<   -1  <<   2  <<   0  <<   1  <<   0 << -0.000745
    <<     1 <<    3  <<  -2  <<   0  <<   1  <<   0 <<  0.000670
    <<     1 <<    2  <<   0  <<  -1  <<  -1  <<   0 << -0.000603
    <<     1 <<    0  <<   0  <<   1  <<  -1  <<   0 << -0.000597
    <<     1 <<   -2  <<   2  <<   1  <<   0  <<   0 <<  0.000542
    <<     1 <<    4  <<  -2  <<  -1  <<   0  <<   0 <<  0.000542
    <<     1 <<   -3  <<   3  <<   0  <<   0  <<  -1 << -0.000541
    <<     1 <<   -2  <<   1  <<   1  <<   0  <<  -1 << -0.000469
    <<     1 <<   -2  <<   3  <<  -1  <<   0  <<  -1 << -0.000440
    <<     1 <<    0  <<  -2  <<   1  <<  -1  <<   0 <<  0.000438
    <<     1 <<   -2  <<  -1  <<   1  <<   0  <<   1 <<  0.000422
    <<     1 <<    4  <<  -2  <<   1  <<   0  <<   0 <<  0.000410
    <<     1 <<   -4  <<   4  <<  -1  <<   0  <<   0 << -0.000374
    <<     1 <<   -4  <<   2  <<   1  <<  -1  <<   0 << -0.000365
    <<     1 <<    5  <<  -2  <<   0  <<   0  <<   0 <<  0.000345
    <<     1 <<    3  <<   0  <<  -2  <<   1  <<   0 <<  0.000335
    <<     1 <<   -5  <<   2  <<   2  <<   0  <<   0 << -0.000321
    <<     1 <<    2  <<   0  <<   1  <<   0  <<   0 << -0.000319
    <<     1 <<    1  <<   3  <<   0  <<   0  <<  -1 <<  0.000307
    <<     1 <<   -2  <<   0  <<   1  <<  -2  <<   0 <<  0.000291
    <<     1 <<    4  <<   0  <<  -1  <<   2  <<   0 <<  0.000290
    <<     1 <<    1  <<  -4  <<   0  <<   0  <<   2 << -0.000289
    <<     1 <<    5  <<   0  <<  -2  <<   0  <<   0 <<  0.000286
    <<     1 <<   -1  <<   0  <<   2  <<   1  <<   0 <<  0.000275
    <<     1 <<   -2  <<   1  <<   0  <<   0  <<   0 <<  0.000271
    <<     1 <<    4  <<  -2  <<   1  <<   1  <<   0 <<  0.000263
    <<     1 <<   -3  <<   4  <<  -2  <<   0  <<   0 << -0.000245
    <<     1 <<   -1  <<   3  <<   0  <<   0  <<  -1 <<  0.000225
    <<     1 <<    3  <<  -3  <<   0  <<   0  <<   1 <<  0.000225
    <<     1 <<    5  <<  -2  <<   0  <<   1  <<   0 <<  0.000221
    <<     1 <<    1  <<   2  <<   0  <<   1  <<   0 << -0.000202
    <<     1 <<    2  <<   0  <<   1  <<   1  <<   0 << -0.000200
    <<     1 <<   -5  <<   4  <<   0  <<   0  <<   0 << -0.000199
    <<     1 <<   -2  <<   0  <<  -1  <<  -2  <<   0 <<  0.000192
    <<     1 <<    5  <<   0  <<  -2  <<   1  <<   0 <<  0.000183
    <<     1 <<    1  <<   2  <<  -2  <<   0  <<   0 <<  0.000183
    <<     1 <<    1  <<  -2  <<   2  <<   0  <<   0 <<  0.000183
    <<     1 <<   -2  <<   2  <<   1  <<   1  <<   0 << -0.000170
    <<     1 <<    0  <<   3  <<  -1  <<   0  <<  -1 <<  0.000169
    <<     1 <<    2  <<  -3  <<   1  <<   0  <<   1 <<  0.000168
    <<     1 <<   -2  <<  -2  <<   3  <<   0  <<   0 <<  0.000162
    <<     1 <<   -1  <<   2  <<  -2  <<   0  <<   0 <<  0.000149
    <<     1 <<   -4  <<   3  <<   1  <<   0  <<  -1 << -0.000147
    <<     1 <<   -4  <<   0  <<   3  <<  -1  <<   0 << -0.000141
    <<     1 <<   -1  <<  -2  <<   2  <<  -1  <<   0 <<  0.000138
    <<     1 <<   -2  <<   0  <<   3  <<   0  <<   0 <<  0.000136
    <<     1 <<    4  <<   0  <<  -3  <<   0  <<   0 <<  0.000136
    <<     1 <<    0  <<   1  <<   1  <<   0  <<  -1 <<  0.000127
    <<     1 <<    2  <<  -1  <<  -1  <<   0  <<   1 <<  0.000127
    <<     1 <<    2  <<  -2  <<   1  <<  -1  <<   0 << -0.000126
    <<     1 <<    0  <<   0  <<  -1  <<  -2  <<   0 << -0.000121
    <<     1 <<    2  <<   0  <<   1  <<   2  <<   0 << -0.000121
    <<     1 <<    2  <<  -2  <<  -1  <<  -1  <<   0 <<  0.000117
    <<     1 <<    0  <<   0  <<   1  <<   2  <<   0 << -0.000116
    <<     1 <<    0  <<   1  <<   0  <<   0  <<   0 << -0.000114
    <<     1 <<    2  <<  -1  <<   0  <<   0  <<   0 << -0.000114
    <<     1 <<    0  <<   2  <<  -1  <<  -1  <<   0 << -0.000114
    <<     1 <<   -1  <<  -2  <<   0  <<  -2  <<   0 <<  0.000114
    <<     1 <<   -3  <<   1  <<   0  <<   0  <<   1 <<  0.000113
    <<     1 <<    3  <<  -2  <<   0  <<  -1  <<   0 <<  0.000109
    <<     1 <<   -1  <<  -1  <<   0  <<  -1  <<   1 <<  0.000108
    <<     1 <<    4  <<  -2  <<  -1  <<   1  <<   0 <<  0.000106
    <<     1 <<    2  <<   1  <<  -1  <<   0  <<  -1 << -0.000106
    <<     1 <<    0  <<  -1  <<   1  <<   0  <<   1 << -0.000106
    <<     1 <<   -2  <<   4  <<  -1  <<   0  <<   0 <<  0.000105
    <<     1 <<    4  <<  -4  <<   1  <<   0  <<   0 <<  0.000104
    <<     1 <<   -3  <<   1  <<   2  <<   0  <<  -1 << -0.000103
    <<     1 <<   -3  <<   3  <<   0  <<  -1  <<  -1 << -0.000100
    <<     1 <<    1  <<   2  <<   0  <<   2  <<   0 << -0.000100
    <<     1 <<    1  <<  -2  <<   0  <<  -2  <<   0 << -0.000100
    <<     1 <<    3  <<   0  <<   0  <<   3  <<   0 <<  0.000099
    <<     1 <<   -1  <<   2  <<   0  <<  -1  <<   0 << -0.000098
    <<     1 <<   -2  <<   1  <<  -1  <<   0  <<   1 <<  0.000093
    <<     1 <<    0  <<  -3  <<   1  <<   0  <<   1 <<  0.000093
    <<     1 <<   -3  <<  -1  <<   2  <<   0  <<   1 <<  0.000090
    <<     1 <<    2  <<   0  <<  -1  <<   2  <<   0 << -0.000088
    <<     1 <<    6  <<  -2  <<  -1  <<   0  <<   0 <<  0.000083
    <<     1 <<    2  <<   2  <<  -1  <<   0  <<   0 << -0.000083
    <<     1 <<   -1  <<   1  <<   0  <<  -1  <<  -1 << -0.000082
    <<     1 <<   -2  <<   3  <<  -1  <<  -1  <<  -1 << -0.000081
    <<     1 <<   -1  <<   0  <<   0  <<   0  <<   2 << -0.000079
    <<     1 <<   -5  <<   0  <<   4  <<   0  <<   0 << -0.000077
    <<     1 <<    1  <<   0  <<   0  <<   0  <<  -2 << -0.000075
    <<     1 <<   -2  <<   1  <<   1  <<  -1  <<  -1 << -0.000075
    <<     1 <<    1  <<  -1  <<   0  <<   1  <<   1 << -0.000075
    <<     1 <<    1  <<   2  <<   0  <<   0  <<  -2 <<  0.000071
    <<     1 <<   -3  <<   1  <<   1  <<   0  <<   0 <<  0.000071
    <<     1 <<   -4  <<   4  <<  -1  <<  -1  <<   0 << -0.000071
    <<     1 <<    1  <<   0  <<  -2  <<  -1  <<   0 <<  0.000068
    <<     1 <<   -2  <<  -1  <<   1  <<  -1  <<   1 <<  0.000068
    <<     1 <<   -3  <<   2  <<   2  <<   0  <<   0 <<  0.000065
    <<     1 <<    5  <<  -2  <<  -2  <<   0  <<   0 <<  0.000065
    <<     1 <<    3  <<  -4  <<   2  <<   0  <<   0 <<  0.000064
    <<     1 <<    1  <<  -2  <<   0  <<   0  <<   2 <<  0.000064
    <<     1 <<   -1  <<   4  <<  -2  <<   0  <<   0 <<  0.000064
    <<     1 <<    2  <<   2  <<  -1  <<   1  <<   0 << -0.000064
    <<     1 <<   -5  <<   2  <<   2  <<  -1  <<   0 << -0.000060
    <<     1 <<    1  <<  -3  <<   0  <<  -1  <<   1 <<  0.000056
    <<     1 <<    1  <<   1  <<   0  <<   1  <<  -1 <<  0.000056
    <<     1 <<    6  <<  -2  <<  -1  <<   1  <<   0 <<  0.000053
    <<     1 <<   -2  <<   2  <<  -1  <<  -2  <<   0 <<  0.000053
    <<     1 <<    4  <<  -2  <<   1  <<   2  <<   0 <<  0.000053
    <<     1 <<   -6  <<   4  <<   1  <<   0  <<   0 << -0.000053
    <<     1 <<    5  <<  -4  <<   0  <<   0  <<   0 <<  0.000053
    <<     1 <<   -3  <<   4  <<   0  <<   0  <<   0 <<  0.000053
    <<     1 <<    1  <<   2  <<  -2  <<   1  <<   0 <<  0.000052
    <<     1 <<   -2  <<   1  <<   0  <<  -1  <<   0 <<  0.000050
    <<     0 <<    2  <<   0  <<   0  <<   0  <<   0 << -0.066607
    <<     0 <<    1  <<   0  <<  -1  <<   0  <<   0 << -0.035184
    <<     0 <<    0  <<   2  <<   0  <<   0  <<   0 << -0.030988
    <<     0 <<    0  <<   0  <<   0  <<   1  <<   0 <<  0.027929
    <<     0 <<    2  <<   0  <<   0  <<   1  <<   0 << -0.027616
    <<     0 <<    3  <<   0  <<  -1  <<   0  <<   0 << -0.012753
    <<     0 <<    1  <<  -2  <<   1  <<   0  <<   0 << -0.006728
    <<     0 <<    2  <<  -2  <<   0  <<   0  <<   0 << -0.005837
    <<     0 <<    3  <<   0  <<  -1  <<   1  <<   0 << -0.005286
    <<     0 <<    0  <<   1  <<   0  <<   0  <<  -1 << -0.004921
    <<     0 <<    2  <<   0  <<  -2  <<   0  <<   0 << -0.002884
    <<     0 <<    2  <<   0  <<   0  <<   2  <<   0 << -0.002583
    <<     0 <<    3  <<  -2  <<   1  <<   0  <<   0 << -0.002422
    <<     0 <<    1  <<   0  <<  -1  <<  -1  <<   0 <<  0.002310
    <<     0 <<    1  <<   0  <<  -1  <<   1  <<   0 <<  0.002283
    <<     0 <<    4  <<  -2  <<   0  <<   0  <<   0 << -0.002037
    <<     0 <<    1  <<   0  <<   1  <<   0  <<   0 <<  0.001883
    <<     0 <<    0  <<   3  <<   0  <<   0  <<  -1 << -0.001811
    <<     0 <<    4  <<   0  <<  -2  <<   0  <<   0 << -0.001687
    <<     0 <<    3  <<  -2  <<   1  <<   1  <<   0 << -0.001004
    <<     0 <<    3  <<  -2  <<  -1  <<   0  <<   0 << -0.000925
    <<     0 <<    4  <<  -2  <<   0  <<   1  <<   0 << -0.000844
    <<     0 <<    0  <<   2  <<   0  <<   1  <<   0 <<  0.000766
    <<     0 <<    1  <<   0  <<   1  <<   1  <<   0 <<  0.000766
    <<     0 <<    4  <<   0  <<  -2  <<   1  <<   0 << -0.000700
    <<     0 <<    3  <<   0  <<  -1  <<   2  <<   0 << -0.000495
    <<     0 <<    5  <<  -2  <<  -1  <<   0  <<   0 << -0.000492
    <<     0 <<    1  <<   2  <<  -1  <<   0  <<   0 <<  0.000491
    <<     0 <<    1  <<  -2  <<   1  <<  -1  <<   0 <<  0.000483
    <<     0 <<    1  <<  -2  <<   1  <<   1  <<   0 <<  0.000437
    <<     0 <<    2  <<  -2  <<   0  <<  -1  <<   0 << -0.000416
    <<     0 <<    2  <<  -3  <<   0  <<   0  <<   1 << -0.000384
    <<     0 <<    2  <<  -2  <<   0  <<   1  <<   0 <<  0.000374
    <<     0 <<    0  <<   2  <<  -2  <<   0  <<   0 << -0.000312
    <<     0 <<    1  <<  -3  <<   1  <<   0  <<   1 << -0.000288
    <<     0 <<    0  <<   0  <<   0  <<   2  <<   0 << -0.000273
    <<     0 <<    0  <<   1  <<   0  <<   0  <<   1 <<  0.000259
    <<     0 <<    1  <<   2  <<  -1  <<   1  <<   0 <<  0.000245
    <<     0 <<    3  <<   0  <<  -3  <<   0  <<   0 << -0.000232
    <<     0 <<    2  <<   1  <<   0  <<   0  <<  -1 <<  0.000229
    <<     0 <<    1  <<  -1  <<  -1  <<   0  <<   1 << -0.000216
    <<     0 <<    1  <<   0  <<   1  <<   2  <<   0 <<  0.000206
    <<     0 <<    5  <<  -2  <<  -1  <<   1  <<   0 << -0.000204
    <<     0 <<    2  <<  -1  <<   0  <<   0  <<   1 << -0.000202
    <<     0 <<    2  <<   2  <<  -2  <<   0  <<   0 <<  0.000200
    <<     0 <<    1  <<  -1  <<   0  <<   0  <<   0 <<  0.000195
    <<     0 <<    5  <<   0  <<  -3  <<   0  <<   0 << -0.000190
    <<     0 <<    2  <<   0  <<  -2  <<   1  <<   0 <<  0.000187
    <<     0 <<    1  <<   1  <<  -1  <<   0  <<  -1 <<  0.000180
    <<     0 <<    3  <<  -4  <<   1  <<   0  <<   0 << -0.000179
    <<     0 <<    0  <<   2  <<   0  <<   2  <<   0 <<  0.000170
    <<     0 <<    2  <<   0  <<  -2  <<  -1  <<   0 <<  0.000153
    <<     0 <<    4  <<  -3  <<   0  <<   0  <<   1 << -0.000137
    <<     0 <<    3  <<  -1  <<  -1  <<   0  <<   1 << -0.000119
    <<     0 <<    0  <<   2  <<   0  <<   0  <<  -2 << -0.000119
    <<     0 <<    3  <<  -3  <<   1  <<   0  <<   1 << -0.000112
    <<     0 <<    2  <<  -4  <<   2  <<   0  <<   0 << -0.000110
    <<     0 <<    4  <<  -2  <<  -2  <<   0  <<   0 << -0.000110
    <<     0 <<    3  <<   1  <<  -1  <<   0  <<  -1 <<  0.000107
    <<     0 <<    5  <<  -4  <<   1  <<   0  <<   0 << -0.000095
    <<     0 <<    3  <<  -2  <<  -1  <<  -1  <<   0 << -0.000095
    <<     0 <<    3  <<  -2  <<   1  <<   2  <<   0 << -0.000091
    <<     0 <<    4  <<  -4  <<   0  <<   0  <<   0 << -0.000090
    <<     0 <<    6  <<  -2  <<  -2  <<   0  <<   0 << -0.000081
    <<     0 <<    5  <<   0  <<  -3  <<   1  <<   0 << -0.000079
    <<     0 <<    4  <<  -2  <<   0  <<   2  <<   0 << -0.000079
    <<     0 <<    2  <<   2  <<  -2  <<   1  <<   0 <<  0.000077
    <<     0 <<    0  <<   4  <<   0  <<   0  <<  -2 << -0.000073
    <<     0 <<    3  <<  -1  <<   0  <<   0  <<   0 <<  0.000069
    <<     0 <<    3  <<  -3  <<  -1  <<   0  <<   1 << -0.000067
    <<     0 <<    4  <<   0  <<  -2  <<   2  <<   0 << -0.000066
    <<     0 <<    1  <<  -2  <<  -1  <<  -1  <<   0 <<  0.000065
    <<     0 <<    2  <<  -1  <<   0  <<   0  <<  -1 <<  0.000064
    <<     0 <<    4  <<  -4  <<   2  <<   0  <<   0 << -0.000062
    <<     0 <<    2  <<   1  <<   0  <<   1  <<  -1 <<  0.000060
    <<     0 <<    3  <<  -2  <<  -1  <<   1  <<   0 <<  0.000059
    <<     0 <<    4  <<  -3  <<   0  <<   1  <<   1 << -0.000056
    <<     0 <<    2  <<   0  <<   0  <<   3  <<   0 <<  0.000055
    <<     0 <<    6  <<  -4  <<   0  <<   0  <<   0 << -0.000051;
    return dnntab;
}

// ocean tide loading
// ---------- 
t_gtriple t_gtide2010::load_ocean(const t_gtime& epoch, const string& site, const t_gtriple& xRec)
{
#ifdef BMUTEX   
   boost::mutex::scoped_lock lock(_mutex);
#endif 
   _mutex.lock();

    Matrix otl, Acjcos, AcjSin;
    Matrix DD, ff, A(2,3), P(2,3); // ff stands for fundamental frequencies
    Matrix it;
    Matrix dnn, dn;
    Matrix votl;
    Matrix fjc(342,1), fjs(342,1), itt(342,6);
    fjc=0.0;
    fjs=0.0;
    double f, fp1, fp2, K, Kp1, Kp2, pp;
	int i1, i2, i1b, i2b;
    double delta = 32.184 + epoch.leapsec();
      
    //                 NOTE !!!
    // pivot indices for each wave in dnn and dnnBLQ - if you modify "dnn" or dnnBLQ (e.g., flipping rows), you HAVE TO change iind accordingly!!!
    // 1. and 2. col determines position of both pivots (if any) in dnn(342,7), 3. and 4. determines position of the pivots (if any) in dnnBLQ (11,6)
    // if iind(i,1)=0 and iind(i,2)=0 and iind(i,1)=iind(i,2) ... pivot from BLQ
    // if iind(i,1)=0 or iind(i,2)=0   						  ... one-side interpolation of the minor wave
    // if iind(i,1)~=0 or iind(i,2)~=0 						  ... two-side interpolation of the minor wave

    int iind[342][4]={
	{0,0,1,1},{0,0,2,2},{0,0,3,3},{0,0,4,4},{4,0,4,0},{3,1,3,1},{3,1,3,1},{3,0,3,0},{1,2,1,2},{1,2,1,2},{3,0,3,0},{3,0,3,0},{1,2,1,2},{3,0,3,0},{4,0,4,0},
	{1,2,1,2},{4,0,4,0},{2,4,2,4},{3,1,3,1},{1,2,1,2},{4,0,4,0},{1,2,1,2},{3,1,3,1},{3,0,3,0},{3,0,3,0},{4,0,4,0},{3,1,3,1},{3,1,3,1},{2,4,2,4},{3,0,3,0},
	{4,0,4,0},{3,0,3,0},{3,1,3,1},{4,0,4,0},{1,2,1,2},{3,0,3,0},{4,0,4,0},{1,2,1,2},{3,0,3,0},{1,2,1,2},{1,2,1,2},{1,2,1,2},{3,1,3,1},{2,4,2,4},{4,0,4,0},
	{3,0,3,0},{3,0,3,0},{4,0,4,0},{3,0,3,0},{3,0,3,0},{1,2,1,2},{3,0,3,0},{4,0,4,0},{4,0,4,0},{1,2,1,2},{1,2,1,2},{1,2,1,2},{3,0,3,0},{3,1,3,1},{4,0,4,0},
	{4,0,4,0},{1,2,1,2},{1,2,1,2},{3,1,3,1},{3,0,3,0},{3,0,3,0},{3,1,3,1},{4,0,4,0},{1,2,1,2},{3,0,3,0},{4,0,4,0},{1,2,1,2},{4,0,4,0},{3,0,3,0},{3,0,3,0},
	{3,0,3,0},{4,0,4,0},{1,2,1,2},{3,0,3,0},{1,2,1,2},{3,1,3,1},{4,0,4,0},{1,2,1,2},{2,4,2,4},{3,0,3,0},{3,0,3,0},{3,0,3,0},{3,0,3,0},{3,0,3,0},{3,1,3,1},
	{3,0,3,0},{4,0,4,0},{1,2,1,2},{3,1,3,1},{3,1,3,1},{1,2,1,2},{3,0,3,0},{3,1,3,1},{3,0,3,0},{3,0,3,0},{4,0,4,0},{3,0,3,0},{3,0,3,0},{3,1,3,1},{3,0,3,0},
	{4,0,4,0},{3,0,3,0},{3,1,3,1},{4,0,4,0},{0,0,5,5},{0,0,6,6},{0,0,7,7},{0,0,8,8},{110,0,5,0},{113,111,8,6},{110,0,5,0},{111,112,6,7},{110,0,5,0},{113,111,8,6},
	{113,0,8,0},{113,0,8,0},{111,112,6,7},{112,110,7,5},{110,0,5,0},{111,112,6,7},{113,0,8,0},{110,0,5,0},{111,112,6,7},{110,0,5,0},{111,112,6,7},{110,0,5,0},
	{110,0,5,0},{111,112,6,7},{110,0,5,0},{112,110,7,5},{110,0,5,0},{113,0,8,0},{111,112,6,7},{113,111,8,6},{110,0,5,0},{111,112,6,7},{113,111,8,6},{110,0,5,0},
	{113,0,8,0},{110,0,5,0},{111,112,6,7},{111,112,6,7},{113,0,8,0},{110,0,5,0},{112,110,7,5},{113,111,8,6},{111,112,6,7},{111,112,6,7},{113,111,8,6},{110,0,5,0},
	{113,0,8,0},{111,112,6,7},{110,0,5,0},{110,0,5,0},{111,112,6,7},{113,111,8,6},{110,0,5,0},{113,0,8,0},{113,111,8,6},{113,111,8,6},{111,112,6,7},{113,0,8,0},
	{110,0,5,0},{113,0,8,0},{113,0,8,0},{110,0,5,0},{110,0,5,0},{113,0,8,0},{110,0,5,0},{110,0,5,0},{113,0,8,0},{110,0,5,0},{111,112,6,7},{110,0,5,0},{111,112,6,7},
	{113,111,8,6},{110,0,5,0},{113,0,8,0},{111,112,6,7},{110,0,5,0},{110,0,5,0},{110,0,5,0},{110,0,5,0},{113,0,8,0},{113,0,8,0},{110,0,5,0},{110,0,5,0},{112,110,7,5},
	{113,111,8,6},{111,112,6,7},{110,0,5,0},{113,0,8,0},{111,112,6,7},{113,0,8,0},{113,0,8,0},{113,111,8,6},{113,111,8,6},{110,0,5,0},{111,112,6,7},{110,0,5,0},
	{110,0,5,0},{111,112,6,7},{110,0,5,0},{110,0,5,0},{111,112,6,7},{111,112,6,7},{110,0,5,0},{111,112,6,7},{113,111,8,6},{113,0,8,0},{110,0,5,0},{113,111,8,6},
	{110,0,5,0},{110,0,5,0},{111,112,6,7},{113,111,8,6},{110,0,5,0},{113,0,8,0},{113,0,8,0},{110,0,5,0},{111,112,6,7},{110,0,5,0},{111,112,6,7},{113,111,8,6},
	{111,112,6,7},{113,0,8,0},{110,0,5,0},{110,0,5,0},{110,0,5,0},{111,112,6,7},{113,111,8,6},{111,112,6,7},{113,0,8,0},{112,110,7,5},{113,111,8,6},{112,110,7,5},
	{110,0,5,0},{113,0,8,0},{113,0,8,0},{112,110,7,5},{113,0,8,0},{113,0,8,0},{110,0,5,0},{110,0,5,0},{112,110,7,5},{111,112,6,7},{110,0,5,0},{113,0,8,0},
	{111,112,6,7},{110,0,5,0},{110,0,5,0},{113,111,8,6},{110,0,5,0},{113,0,8,0},{110,0,5,0},{113,0,8,0},{110,0,5,0},{113,111,8,6},{0,0,9,9},{0,0,10,10},{0,0,11,11},
	{266,0,11,0},{264,0,9,0},{264,0,9,0},{266,265,11,10},{265,264,10,9},{264,0,9,0},{266,0,11,0},{265,264,10,9},{264,0,9,0},{264,0,9,0},{266,265,11,10},
	{265,264,10,9},{264,0,9,0},{265,264,10,9},{266,265,11,10},{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0},{266,265,11,10},{265,264,10,9},{264,0,9,0},{264,0,9,0},
	{264,0,9,0},{265,264,10,9},{266,265,11,10},{266,265,11,10},{265,264,10,9},{265,264,10,9},{265,264,10,9},{266,0,11,0},{266,265,11,10},{266,0,11,0},{266,0,11,0},
	{265,264,10,9},{264,0,9,0},{264,0,9,0},{266,265,11,10},{265,264,10,9},{264,0,9,0},{265,264,10,9},{264,0,9,0},{266,265,11,10},{264,0,9,0},{265,264,10,9},
	{265,264,10,9},{264,0,9,0},{266,265,11,10},{265,264,10,9},{264,0,9,0},{264,0,9,0},{266,0,11,0},{264,0,9,0},{265,264,10,9},{264,0,9,0},{264,0,9,0},{264,0,9,0},
	{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0},{266,265,11,10},{264,0,9,0},{264,0,9,0},{264,0,9,0},{266,265,11,10},
	{265,264,10,9},{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0},{264,0,9,0}
	};
    
	Matrix dnnBLQ(11,6);
	dnnBLQ	<< 2 <<  0 <<  0 <<  0 << 0 << 0
			<< 2 <<  2 << -2 <<  0 << 0 << 0
			<< 2 << -1 <<  0 <<  1 << 0 << 0
			<< 2 <<  2 <<  0 <<  0 << 0 << 0
			<< 1 <<  1 <<  0 <<  0 << 0 << 0
			<< 1 << -1 <<  0 <<  0 << 0 << 0
			<< 1 <<  1 << -2 <<  0 << 0 << 0
			<< 1 << -2 <<  0 <<  1 << 0 << 0
			<< 0 <<  2 <<  0 <<  0 << 0 << 0
			<< 0 <<  1 <<  0 << -1 << 0 << 0
			<< 0 <<  0 <<  2 <<  0 << 0 << 0 ;
      
    // load dMJD
    double mjd_i=epoch.dmjd(); // [day]
    
    // input for the atronomical argument for each epoch
    DD=DD_arguments(mjd_i, delta);
      
   // load all waves for the interpolation
    dnn=dnn_table();    
    
    // load BLQ for a particular site
	t_gtriple dxyz(0.0,0.0,0.0);  
	if( !_gotl || _gotl->data(otl, site) < 0 ){
    if(_log) _log->comment(0, "gtide2010", "WARNING: Site " + site + " not found in ocean tide BLQ file!");
    _mutex.unlock(); return dxyz;
  }

	otl.SubMatrix(4,6,1,11)=otl.SubMatrix(4,6,1,11)*D2R;
	
	// ff is the same for all waves and all epochs
	ff = DD.submatrix(7,12,1,1);
	
	Kp2 = 0.0;
	Kp1 = 0.0;
	pp = 0.0;
	
    for (int i = 0; i <= 341; i++ ){
		
		dn = dnn.SubMatrix(i+1,i+1,1,6);
		K = dnn(i+1,7);
		f = DotProduct(dn,ff);
		i1 = iind[i][0]; // indeces wrt dnn, +1 because dnn and otl is Newmat Matrix
		i2 = iind[i][1];
		i1b = iind[i][2]; // indeces wrt BLQ
		i2b = iind[i][3];
		A = 0.0;
		P = 0.0;
			
		// BLQ - pivots not needed at all
		if (i1==0 && i2==0 ){
			fp1 = 0.0; // frequency of the first pivot
			fp2 = 0.0; // frequency of the seoond pivot
			Kp1 = K;
			Kp2 = 1.0; // because of denominator
			pp = 0.0;
				
			A(1,1) = otl(1,i1b);
			A(1,2) = otl(2,i1b);
			A(1,3) = otl(3,i1b);
			
			P(1,1) = otl(4,i1b);
			P(1,2) = otl(5,i1b);
			P(1,3) = otl(6,i1b); //% up,ew,ns    		
		}
		
		// One-side interpolation ... pivot 2 non-existing
		if (i1!=0 && i2==0){
 
 		    fp1 = DotProduct(dnn.SubMatrix(i1,i1,1,6),ff);
			fp2 = 0.0; // frequency of the second pivot
			Kp1 = dnn(i1,7);
			Kp2 = 1.0;  // because of denominator
			pp = 0.0;
			
			A(1,1) = otl(1,i1b);
			A(1,2) = otl(2,i1b);
			A(1,3) = otl(3,i1b);
			
			P(1,1) = otl(4,i1b);
			P(1,2) = otl(5,i1b);
			P(1,3) = otl(6,i1b); //% up,ew,ns   		 
		}

		// One-side interpolation ... pivot 1 non-existing
		if (i1==0 && i2!=0){
 			fp2 = 0.0; // frequency of the first pivot
 		    fp1 = DotProduct(dnn.SubMatrix(i2,i2,1,6),ff);
			Kp2 = 1.0;  // because of denominator
			Kp1 = dnn(i2,7);
			pp = 0.0;
			
			A(1,1) = otl(1,i2b);
			A(1,2) = otl(2,i2b);
			A(1,3) = otl(3,i2b);
			
			P(1,1) = otl(4,i2b);
			P(1,2) = otl(5,i2b);
			P(1,3) = otl(6,i2b); //% up,ew,ns   		 
		}		
		
		// Two-side interpolation
		if (i1!=0 && i2!=0 ){
			fp1 = DotProduct(dnn.SubMatrix(i1,i1,1,6),ff); // frequency of the first pivot
			fp2 = DotProduct(dnn.SubMatrix(i2,i2,1,6),ff); // frequency of the second pivot
			
			// flip all if fp1 > fp2
			if (fp1 > fp2){
				i1 = iind[i][1]; // indeces wrt dnn, +1 because dnn and otl is Newmat Matrix
				i2 = iind[i][0];
				i1b = iind[i][3]; // indeces wrt BLQ
				i2b = iind[i][2];
				fp1 = DotProduct(dnn.SubMatrix(i2,i2,1,6),ff); // frequency of the first pivot
			    fp2 = DotProduct(dnn.SubMatrix(i1,i1,1,6),ff); // frequency of the second pivot
			}	
			Kp1 = dnn(i1,7); // amplitude of the first pivot
		    Kp2 = dnn(i2,7); // amplitude of the second pivot	
		    pp = (f - fp1)/(fp2 - fp1);
		    
		    A(1,1) = otl(1,i1b); P(1,1) = otl(4,i1b);
			A(1,2) = otl(2,i1b); P(1,2) = otl(5,i1b);
			A(1,3) = otl(3,i1b); P(1,3) = otl(6,i1b); //% up,ew,ns
			
			A(2,1) = otl(1,i2b); P(2,1) = otl(4,i2b);
			A(2,2) = otl(2,i2b); P(2,2) = otl(5,i2b);
			A(2,3) = otl(3,i2b); P(2,3) = otl(6,i2b); //% up,ew,ns
		}
	
		// ==== Conditions for Phase and Amplitude
	
		 //if (K>0.0 && dn(1,1)==0){ // ADMINTF degree 0 ... likely not needed
            // P = P + G_PI;       }
         		 
         if (K>0.0 && abs(dn(1,1)-1.0)<0.001){ //ADMINTF degree 1
            P = P + G_PI/2.0;
            A = -1.0*A;          } 
            
         if (K<0.0 && abs(dn(1,1)-1.0)<0.001){ //ADMINTF degree 1
            P = P - G_PI/2.0;
            A = -1.0*A;          }             

        // table itt - interpolation of all tides (for major tides - pivots K/Kp1=1 and pp=0)
        for(unsigned ii = 1; ii <= 3; ii++){
			itt(i+1,ii) = A(1,ii)*cos(P(1,ii))*K/Kp1*(1.0-pp) + A(2,ii)*cos(P(2,ii))*K/Kp2*pp;
			itt(i+1,ii+3)=A(1,ii)*sin(P(1,ii))*K/Kp1*(1.0-pp) + A(2,ii)*sin(P(2,ii))*K/Kp2*pp;		
		}	
        	
    } // ends loop over tide waves
     	
	for(unsigned j = 1; j <= 342; j++){
	 	fjc(j,1) = cos(null2pi_rad(DotProduct(dnn.submatrix(j,j,1,6), DD.submatrix(1,6,1,1))));
		fjs(j,1) = sin(null2pi_rad(DotProduct(dnn.submatrix(j,j,1,6), DD.submatrix(1,6,1,1))));
      } 
	   
	 double up, ns, ew;
     up = DotProduct(fjc,itt.Column(1)) + DotProduct(fjs,itt.Column(4));
     ns = DotProduct(fjc,itt.Column(2)) + DotProduct(fjs,itt.Column(5));
     ew = DotProduct(fjc,itt.Column(3)) + DotProduct(fjs,itt.Column(6));

//	cout << setw(12) <<up <<" "<< ns << " " << ew << endl; 
	
   t_gtriple dneu;
   dneu[0] = -ns;
   dneu[1] = -ew;
   dneu[2] = up;    

   t_gtriple ell;
   xyz2ell(xRec, ell, false);
   neu2xyz(ell, dneu, dxyz);
   
   _mutex.unlock(); return dxyz;
}


// atmospheric tide loading
// ----------
t_gtriple t_gtide2010::load_atmosph()
{
#ifdef BMUTEX   
   boost::mutex::scoped_lock lock(_mutex);
#endif 
   _mutex.lock();

  t_gtriple dxyz(0.0,0.0,0.0);

  _mutex.unlock(); return dxyz;
}

// PROTECTED
//======================   
} // namespace
 
